基因组比对绘图 -- 8703项目

一.单拷贝基因共线性图(单条染色体)

1.准备目标区间的pep序列
## 对每个样本进行提取
CHROM=`awk -F ">"  'NR==1 {print $2}'  chromosome_4A.fasta`  && START=`awk 'NR==1 {print $1}' 12m.bed`  &&  END=`awk 'NR==1 {print $2}' 12m.bed`    && awk -v chrom="$CHROM" -v start="$START" -v end="$END" '{    if ($1 == chrom && $4 >= start && $5 <= end) {        print;    }}' chromosome_4A.gff   > 12m.gff && gffread -x cds.fa -g chromosome_4A.fasta 12m.gff  &&  ir.py -i cds.fa -o pep.fa  -s7   -n 1
## 复制备用
for i in *;do cp $i/pep.fa   gene_plot/$i.pep.fa;done

2.运行聚类
perl /share/nas6/xul/program/fungi/tree/EasyTree.pl -id pep # pep是目录，里面放各自的蛋白序列
## 获得Orthogroups_SingleCopyOrthologues.txt、Orthogroups.txt，路径TREE/format_seq/OrthoFinder/Results_Aug13/Orthogroups

3.绘图
for i in *;do cp $i/12m.gff   gene_plot/TREE/format_seq/OrthoFinder/Results_*/Orthogroups/$i.12m.gff;done
cat Orthogroups_SingleCopyOrthologues.txt  | while read i ;do grep $i Orthogroups.txt  ;done | cut -d ":" -f 2 |sed 's/^ //g' > tmp.txt

## 调整画图的列数和顺序，注意这里索引从0开始
perl -lane 'print join" ",@F[4,5,6,1,2,3]' tmp.txt  | less  > tmp2.txt # 按第5,6,7,2,3,4列重排

## 可以不调整画图的列数，不使用tmp2
perl /share/nas6/xul/program/nuclear/other/get_single_copy_gene_pos.pl -i tmp2.txt  -g  *.12m.gff*  && realpath tmp_collinearity/input.cfg
# 这一步报错可能是因为tmp文件每个基因多了transcript标签，要删掉

## 6种颜色
perl /share/nas6/xul/program/fungi/colinearity/draw_collinearity.pl -i tmp_collinearity/input.cfg  -l tmp_collinearity/*.txt -o synteny.svg && svg2xxx synteny.svg ./ -t png -dpi 300
## 15种颜色，修改线宽为5
perl /share/nas1/yuj/program/fungi/colinearity/draw_collinearity.pl -i tmp_collinearity/input.cfg  -l tmp_collinearity/*.txt2 -o synteny.svg && svg2xxx synteny.svg ./ -t png -dpi 300

通过全部减去720000000来完成部分绘图

二.基因组SV可视化绘图

1.0 区间寻找

1.
## 中国春
cp_add.py -i chromosome_*A.fasta -p 740000001-754227511:+ -sn 4A_12M.fa  > log ; sed -i '1 i\>4A_12M'  *4A_12M.fa && rm log

2.绘图查找
for i in {aegilopoides_TA299,dicoccoides_WEW_v2.1,durum,monococcum_TA10622,urartu_G1812} ;do cd $i && nucmer  -t 15    -l 300 -c 1000 --mum     ../4A_12M.fa  chromosome_*A.fasta  -p $i && delta-filter -m -i 85 -l 15000 -o 85 $i.delta -1>$i.filtered.delta && show-snps -Clr $i.filtered.delta > $i.filtered.delta.snps && mummerplot -p $i $i.filtered.delta -t postscript  && ps2pdf $i.ps $i.pdf && convert -density 300 $i.pdf $i.png & done

3.根据比对结果填写取出的bed文件
cat */*bed  && echo -e "\t$(awk 'NR==1 {print $3}' */*bed)" > 12m.bed && rm *4A_12M.fa  ; cpath=`pwd`&& awk 'NR==4 {print $0}' ../plot/`basename  $cpath`.filtered.delta

4.提取序列
start=`awk   'NR==1 {print $1}' 12m.bed`   && end=`awk   'NR==1 {print $2}' 12m.bed` &&  cp_add.py -i chromosome_4A.fasta -p  $start"-"$end":+"  -sn 4A_12M.fa  > log ; sed -i '1 i\>4A_12M'  *4A_12M.fa  && rm log && ir.py -l -i 4A_12M.fa && ll

1.1提取

1.查看长度
seqkit fx2tab --length --name --header-line  *.fna > size # id（包括染色体）<TAB>长度

2.提取某染色体序列和gff
fa=*.f*a* && gff=*.gff* && seqkit fx2tab --length --name --header-line  $fa > size 
chr=4A && grep $chr size | awk '{print  $1"\t"0"\t"$NF}'  > $chr.bed && bedtools getfasta -fi $fa -bed $chr.bed -fo $chr.fasta && tail -n 1 $chr.fasta > noid.$chr.fasta && fold -w 80 noid.$chr.fasta > chromosome_${chr}.fasta && sed -i '1 i\>4A' chromosome_4A.fasta && grep `cat $chr.bed|awk '{print  $1}'` $gff  > chromosome_${chr}.gff && mv chromosome_* ../

3.提取染色体序列末尾
length=`cut -f 3 $chr.bed` && cd ../ && cp_add.py -i chromosome_4A.fasta -p 740000001-${length}:+ -sn 4A_12M.fa  > log ; sed -i '1 i\>4A_12M'  *4A_12M.fa && rm log

1.2比对画图

1.2.1 流程minimap2

/share/nas6/xul/program/nuclear/syri/multi-genome_syri_plot.pl

perl /share/nas1/yuj/program/nuclear/syri/multi-genome_syri_plot.pl -i ref_4A_12M.fa dicoccoides_4A_12M.fa xxx.fa  aaa.fa

## 分步进行
# Using minimap2 for generating alignment. Any other whole genome alignment tool can also be used.
1.全基因组序列比对
minimap2 -ax asm5 --eqx refgenome qrygenome > out.sam

2.基于比对结果检测结构变异
python3 $PATH_TO_SYRI -c out.sam -r refgenome -q qrygenome -k -F S
# or
samtools view -b out.sam > out.bam
python3 $PATH_TO_SYRI -c out.bam -r refgenome -q qrygenome -k -F B

3.结构变异检测结果绘图
python3 $PATH_TO_PLOTSR syri.out refgenome qrygenome -H 8 -W 5

1.2.2 mummer

1.基因组比对
# nucmer  -t 15 --maxmatch -c 100 -b 500 -l 50 refgenome qrygenome 
# Whole genome alignment. Any other alignment can also be used. -maxmatch 
# 会识别所有比对，搭配“-m -i 85 -l 15000 -o 85 out.delta -1”绘图会比较好看
nucmer  -t 15    -l 300 -c 1000 --mum refgenome qrygenome # 更快

2.比对结果过滤
delta-filter -m -i 85 -l 15000 -o 85 out.delta -1> out.filtered.delta # Remove small and lower quality alignments 筛选长度>15000bp，匹配度>85% 的匹配情况
-m 会删除冗余比对。
-i 指定最小的alignment相似性阈值
-l 注意，这里是字母小写的L，指定最小的alignment长度
-o 和-r,-q相关，可以理解为alignment coverage
-1 注意，这里是数字1，指定是否进行一对一的比对，一个位置(subject或query上)只找一个最佳的比对。特别是对大的基因组一定要加这个选项，否则会异常慢

show-snps -Clr out.filtered.delta > out.filtered.delta.snps
mummerplot -p out out.filtered.delta -t postscript # ps格式图片
ps2pdf out.ps out.pdf && convert -density 300 out.pdf out.png

for i in *.fa;do echo ${i%_4A_12M.fa};nucmer  -t 15    -l 300 -c 1000 --mum  ../CS2.1_4A_12M.fa $i -p ${i%_4A_12M.fa}   && delta-filter -m -i 85 -l 15000 -o 85   ${i%_4A_12M.fa}.delta -1>  ${i%_4A_12M.fa}.filtered.delta && show-snps -Clr ${i%_4A_12M.fa}.filtered.delta  > ${i%_4A_12M.fa}.filtered.delta.snps && mummerplot -p ${i%_4A_12M.fa}  ${i%_4A_12M.fa}.filtered.delta  -t postscript  && ps2pdf ${i%_4A_12M.fa}.ps  ${i%_4A_12M.fa}.pdf  && convert -density 300   ${i%_4A_12M.fa}.pdf  ${i%_4A_12M.fa}.png &  done

3.获得每个 alignment 的位置
show-coords -THrd out.filtered.delta > out.filtered.coords      # Convert alignment information to a .TSV format as required by SyRI

4.结构变异检测
syri -c out.filtered.coords -d out.filtered.delta -r refgenome -q qrygenome --prefix 4A_12M --nc 8 # -nc 核心数，默认1

5.结果绘图
plotsr 4A_12Msyri.out refgenome qrygenome  -H 8 -W 12 -o 4A_12M.pdf # 简简单单
/share/nas1/yuj/software/micromamba/envs/py3.8/bin/plotsr --genomes genomes.txt -o output_plot.pdf -H 8 -W 12  --sr ref_4A_12M_dicoccoides_4A_12Msyri.out # --genomes 是绘图样式参数，里面有ref qry路径

结果文件

syri.out #检测序列变异的文本信息
syri.vcf #检测序列vcf
syri.summary #变异类型统计文件
syri.pdf #输出图

syri.out
顺序依次为：参考染色体/起始位置/结束位点/参考位点/等位位点/比对染色体(query)/起始位点(query)/结束位点(query)/类型+顺序(ref)/类型+顺序(query)/注释类型/复制状态

SYRI https://www.jianshu.com/p/9eed9bb21cd9

输入文件：
-c INFILE：包含比对坐标的文件。这个文件可能是从一个比对软件（比如MUMmer）生成的，其中包含了两个基因组之间的比对信息，包括每个比对的起始位置、终止位置等信息。
-r REF：作为比对参考的基因组A的文件。
-q QRY：作为比对查询的基因组B的文件。
-d DELTA：MUMmer生成的.delta文件。

可选参数：
-F {T,S,B}：指定输入文件的类型。可以是表格文件（T），SAM文件（S）或BAM文件（B）。默认为表格文件。
-k：保留中间输出文件。默认为False。
--log {DEBUG,INFO,WARN}：设置日志级别。可以选择DEBUG（调试）、INFO（信息）、WARN（警告）等级别。默认为INFO。
--lf LOG_FIN：指定日志文件的名称。默认为"syri.log"。
--dir DIR：指定工作目录的路径。默认为当前目录。
--prefix PREFIX：指定输出文件名前缀。默认为空。
--seed SEED：用于生成随机数的种子。默认为1。
--nc NCORES：指定并行计算时使用的核心数量。最大值是染色体的数量。默认为1。
--novcf：不将所有文件合并为一个输出文件。默认为False。
-f：过滤掉低质量的比对。默认为True。

结构变异检测：
--nosv：设置为True以跳过结构变异的识别。默认为False。
--nosnp：设置为True以跳过在比对中识别SNP/Indel。默认为False。
--all：设置为True以使用重复区域（duplications）进行变异识别。默认为False。
--allow-offset OFFSET：允许碱基对（base pairs）重叠的数量。默认为5，表示如果两个变异区域之间的重叠小于等于5个碱基对，则允许识别为同一变异。
--cigar：设置为True以使用CIGAR字符串来查找SNP/Indel。默认为False，表示不使用CIGAR字符串。通常在使用除了nucmers之外的比对器生成的比对结果时，需要设置为True。
-s SSPATH：指定show-snps工具的路径，show-snps是MUMmer软件包中的一个工具，用于从.delta文件中提取SNP信息。默认为show-snps。

结构重排识别：
--nosr：设置为True以跳过结构重排事件的识别。默认为False，表示执行结构重排事件的识别。
--tdgaplen TDGL：多重比对转座或重复（TD）的两个比对之间允许的最大间隙长度。较大的值会增加TD识别的灵敏度，但也会增加运行时间。默认为500000。
-b BRUTERUNTIME：限制Brute Force方法运行时间的阈值（以秒为单位）。较小的值会使算法运行更快，但可能会对准确性产生边际影响。在一般情况下，可能不需要设置这个参数。默认为60。
--unic TRANSUNICOUNT：选择转座时所需的唯一碱基对数。较小的值会更好地选择较小的转座，但可能会增加时间并降低准确性。默认为1000。
--unip TRANSUNIPERCENT：选择转座时所需的唯一区域百分比。值应在(0,1]范围内。较小的值会选择与其他区域更重叠的转座。默认为0.5。
--inc INCREASEBY：添加另一个比对到转座簇解决方案所需的最小分数增加。默认为1000。
--no-chrmatch：如果两个基因组的染色体ID不相等，则设置为True，不允许SyRI自动匹配染色体ID。默认为False。

mummer

show-snps 用于显示两样本的snp信息
show-aligns 用于显示比对，可以单独列出每个序列的比对情况。
show-coords 用于显示比对坐标，用于后续共线性分析定制化绘图
show-diff显示大的染色体变化 倍增 重排或者直接使用dnadiff软件一步生成，结果非常详细，还有一个report报告文件
dnadiff可以直接加-d接delta格式的结果（/opt/software/mummer-3.9.4alpha/dnadiff -d \<outpfix\>.delta），或者更方便直接接两条序列即可，非常方便好用。

plotsr

python /path/to/plotsr syri.out /path/to/refgenome /path/to/qrygenome

positional arguments:
  reg                   syri.out file generated by SyRI
  r                     path to reference genome
  q                     path to query genome

optional arguments:
  -h, --help            show this help message and exit
  -s S                  minimum size of a SR to be plotted
  -R                    Create ribbons
  -f F                  font size
  -H H                  height of the plot
  -W W                  width of the plot
  -o {pdf,png,svg}      output file format (pdf, png, svg)
  -d D                  DPI for the final image
  -b {agg,cairo,pdf,pgf,ps,svg,template}
                        Matplotlib backend to use

1.3其他操作

1.seqkit
seqkit stat  *.fna # 统计总长
seqkit seq -rp ly.chr1.fasta > ly.chr1-rp.fasta # 反向互补

2.统计长度
samtools faidx *.fna && cat *.fai | cut -f 1,2 > size # id（只有序列号）<TAB>长度

3.提取某染色体 # https://www.jianshu.com/p/e834562ef52c/
$ vim 4A.bed #每条Scaffold的长度信息，一般基因组注释文件gff含有这些信息
NC_057803.1	0	754227511
NC_057803.1	0	754227511
$ bedtools getfasta -fi *.fna -bed 4A.bed -fo 4A.fasta
grep "^>" -v 4A.fasta | awk '{ ORS = ""; $1 = $1; print $0}' > chrUn.fasta # 去掉名字
fold -w 80 chrUn.fasta > chromosome_4A.fasta # 控制行宽80输出
sed -i '1 i\>NC_057803.1' chromosome_4A.fasta # 在第一行加">ID"

4.切分为两段 # https://www.jianshu.com/p/e834562ef52c/
$ vim chr4A_2parts.bed # 这里每条染色体从哪里断开是随意设置的，有一点点不严谨
NC_057803.1	0	400000000
NC_057803.1	400000000	754227511
$ bedtools getfasta -fi *.fna -bed chr4A_2parts.bed -fo chr4A_2parts.fasta
grep "^>" chr4A_2parts.fasta
fold -w 80 chr4A_2parts.fasta > chromosome_4A_2parts.fasta

三.xul

bash /share/nas2/yuj/project/2024/genome/GP-20240704-8703_20240717/region.sh  chromosome_4A.gff chr4A  740000000  800000000  >  chr4AL.gff

个性分析 2DL 特有序列
        1.获取CDS序列
                gffread  -x yzm_chr2AL.fa -g GCA_900231445.1_Svevo.v1_genomic.fna  chr2AL.gff

        2.比对
                makeblastdb  -dbtype nucl -in final.cds.fa
                for i in chym_chr2E.cds.fa dm_chr2H.cds.fa  hm_2R.cds.fa ;do blastn  -db final.cds.fa -query $i  -outfmt 5 > $i.blast.xml ;done
                for i in *.blast.xml ;do perl /share/nas6/pub/pipline/genome-assembly-seq/mitochondrial-genome-seq/current/annotation/src/blast_parser.pl -tophit 1 -topmatch 1 -eval 1e-5 -m 7 $i > `basename $i .blast.xml`.blast.tophit.xls;done 

        3.整理结果
                perl /share/nas6/xul/program/nuclear/other/get_homo_cds_from_blastn_tophit.pl -i *.tophit.xls -g *.gff -f final.cds.fa -o jjm_and_yzm_L.homo.pos.xls  > jjm_and_yzxm_L.homo.xls
                #perl /share/nas6/xul/program/nuclear/other/get_homo_cds_from_blastn_tophit.pl -i *.tophit.xls  ../chr2DL.cds.filter.fa.blast.tophit.xls -g *.gff3  ../chr2DL.gff -f ../2VL.cds.filter.fa   -o  homo.pos.xls

                venn：
                        head -7 ptxm.homo.xls  | sed 's/#//' > ptxm.homo.venn
                        perl /share/nas6/xul/program/nuclear/gene_family/veen_draw.pl  -i  jjm_and_yzxm.homo.venn -o jjm_and_yzxm.homo.venn.svg

        4.共线性图
                perl /share/nas6/xul/program/nuclear/other/get_link_from_homo.pos.pl -i all.homo.pos.xls 

                cp /share/nas6/xul/program/nuclear/other/ticks.conf /share/nas6/xul/program/nuclear/other/ideogram.conf .
                
                /share/nas6/xul/soft/circos/circos-0.69-6/bin/circos  -conf circos.conf && svg2xxx -t pdf circos.svg